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Abstract 

An equivalent domain integral (EDI) method for calculating J — integrals for 
two-dimensional cracked elastic bodies is presented. The details of the method 
and its implementation axe presented for isoparametric elements. The total and 
product integrals consist of the sum of an area or domain integral and line integrals 
on the crack faces. The line integrals vanish only when the crack faces are traction 
free and the loading is either pure mode I or pure mode II or a combination of both 
with only the square-root singular term in the stress field. The EDI method gave 
accurate values of the J— integrals for two mode I and two mixed mode problems. 
Numerical studies showed that domains consisting of one layer of elements axe 
sufficient to obtain accurate J— integral values. Two procedures for separating the 
individual modes from the domain integrals are presented. The procedure that 
uses the symmetric and antisymmetric components of the stress and displacement 
fields to aalculate the individual modes gave accurate values of the integrals for all 
problems analyzed. The EDI method when applied to a problem of an interface 
crack in two different materials showed that the mode I and mode II components 
are domain dependent while the total integral is not. This behavior is caused by 
the presence of the oscillatory part of the singularity in bimaterial crack problems. 
The EDI method, thus, shows behavior similar to the virtual crack closure method 
for bimaterial problems. 



INTRODUCTION 


Stress-intensity factors or strain-energy release rates are used to characterize 
the severity of the stress and strain fields at cracks or delaminations in elastic 
bodies. Numerous procedures exist in the literature to compute the stress-intensity 
factors or strain-energy release rates by analytical and numerical methods. For 
cracks in elastic-plastic materials, the J — integral was proposed [1-3]. As originally 
formulated, the J— integral is a closed contour integral of the strain-energy density 
and work done by tractions around the crack tip. The J— integral is a path- 
independent parameter and is equivalent to the rate of change of total potential 
energy with reference to the crack length. For elastic bodies the J— integral is 
equivalent to the strain-energy release rate. 

Although the J— integral is a potential parameter to define the severity of the 
crack tip in elastic and elastic-plastic materials, it is cumbersome to compute 
in finite element analyses. Therefore, alternate forms of computing this integral 
have been proposed [4-11]. In references 6 and 7, procedures based on the virtual 
crack extension technique [4,5] were proposed to compute the strain-energy release 
rates (and hence the J— integral). In references 8 through 10 line J— integrals 
were converted to equivalent area or domain integrals, hence the name equivalent 
domain integral (EDI). The conversion of line integrals to domain integrals is 
very attractive because all the quantities necessary for computation of the domain 
integrals are readily available in a finite element analysis. < 

The potential of the EDI method appears to be in the calculation of the 
J— integrals (and the strain-energy release rates, G, for elastic materials) for three- 
dimensional (3D) crack configurations. With this method the J— or G- values 
can be computed at any point on the crack front in a straightforward manner. 
The force method [12] and virtual crack closure method [13], used for 3D crack 
configurations, require that the finite element mesh be normal to the crack front. 
In contrast, the EDI method does not require such customized modeling. The EDI 
method also appears to be a powerful tool to calculate the J— integrals for elastic- 
plastic materials and problems involving other material nonlinearities. Although 
the long term objectives of the present research are to apply the EDI method to 3D 
elastic and elastic-plastic problems, issues such as the mode separation in mixed 
mode situations and problems where the crack faces are subjected to external 
loading need to addressed. These issues are studied for elastic two-dimensional 
crack configurations. Therefore, the purpose of this report is to address these 
two specific issues. First, a detailed formulation of the EDI method applicable to 
mixed mode problems under general loading is presented. Second methods for the 
separation of modes in mixed mode crack problems are studied. 
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As this report is intended to explore the method, the EDI method and 
its numerical implementation in the context of a finite element analysis with 
isoparametric elements is also presented. The method is then applied to two mode 
I problems and two mixed mode problems to study the sensitivity of fracture 
parameters to variables used in the EDI algorithm and to evaluate the methods 
for the separation of modes in mixed mode problems. 

SYMBOLS 


b 

A 

E 

E\, Ei 

JhJii 

J AS x\ > J AS 
JSx^JSxi 

Jxi , ■IlJ 

[J] 

N i 

n 

Uj 

P 

S 

« 1,«2 

U IAS, U 2AS 


half-length or length of crack 
half-width or width of plate 

domain considered in the domain integral evaluation 
Young’s modulus of the homogeneous material 
Young’s moduli of materials 1 and 2 in a bimaterial plate 
mode I and mode II J— integrals 
Antisymmetric J and product integrals 
Symmetric J and product integrals 
J and product integrals 
Jacobian matrix of an isoparametric element 
shape function for the j** 1 node of an isoparametric element 
normal vector to the contour T 
component of the normal n in the j th direction 

remote uniform applied stress or uniform crack face pressure loading 
S-function defined in Eq. (7) 

displacements in the x\— and a^-directions, respectively 

antisymmetric components of displacements in the x\ — and ^-directions, 
respectively 
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«1S>«2S 

W 


Xl,X2 


r,r 0 ,r! 


‘a 

V 


V\,V 2 


symmetric components of displacements in the xi— and X2-directions, 
respectively 

strain-energy density 

Cartesian coordinates 

contours around the crack tip 

strain tensor 

Poisson’s ratio of a homogeneous material 

Poisson’s ratios of materials 1 and 2 of a bimaterial plate 

stress tensor 


Superscripts: 

T denotes the transpose of a matrix or column vector 


CONTOUR AND DOMAIN INTEGRALS 

Consider a crack in a plate subjected to an arbitrary remote loading with an 
arbitrary closed contour T around the crack tip, as shown in Figure 1. Then the 
J — integral is defined in the absence of any body forces as [1-3] 


J Xle — 


L 


Q dr 


( 1 ) 


where k = 1,2 and 


Q — [Writ °ij n j\ 

In Eq. (2) W is the total strain-energy density defined as 


W = 





( 2 ) 


( 3 ) 
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In Eqs. (2) and (3) er,y is the stress tensor, e,j is the total strain (which is the sum 
of the elastic and plastic strains) and nj is the component in the j th direction of 
the vector n normal to the contour I\ The indices t ,j, and k take the values 1 
and 2 for two-dimensional (2D) problems. The x\ and £2 represent the directions 
along and normal to the crack line, respectively. 

The integrals J Xl and J X2 are two path-independent integrals that define the 
total amount of the energy flux leaving the contour T in the two directions x\ 
and £ 2 , respectively. J Xl is usually termed the J — integral and J X2 is called the 
product integral. 

As previously mentioned, the line integrals of Eq. (1) are cumbersome to 
compute in a finite element analysis. In contrast, an area or domain integral is very 
convenient to compute in a finite element analysis. Therefore, an alternate form 
of Eq. (1), called the equivalent domain integral (EDI), was proposed [8,10,11]. A 
more detailed formulation of the EDI applicable for mixed mode problems and its 
numerical implementation is presented here. 

Consider two contours To(OABCO) and Y\{ODEFO) around the crack tip 
as shown in Figure 1(b). The two contours will enclose an area DEFCBAD. By 
multiplying the integral over To by unity and the integral over Ti by zero, J Xk 
can be expressed as 


'*k 


= l/ QdT - of Q dT 

J r 0 JTi 


(4) 


This manipulation is performed to convert the line integrals into an area or a 
domain integral. Equation (4) can be expanded as 


' Xk 


Q dT 


= 1 f Q dT - 0 [ 

Jabc Jdef 

+ 1 f QdT + 1 f Q dT 
Jco JOA 


(5) 


An arbitrary but continuous function *9 (£ 1 , £ 2 ) is introduced that has the 
property 


^(*l> x 2) — 1 
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on r ABC an< i 


S’(xi,x 2 ) = 0 


on Tqef. Using this S-function, Eq. (5) can be written as 


T x . = f Q S dT - f Q S dT + f QdT + [ Q dT 

J ABC Jdef Jco Joa 


= f - [ QSdT - f QSdrl + f QdT + / QdT 

1 JCBA JDEF J JCO Joa 


A closed contour integral along DEFCBAD can be defined easily by adding 
and subtracting the line integrals on the crack faces FC and AD as 


-L 


- f Q S dT - f Q S dT - f Q S dT 
J cba Jdef Jfc 

+ f Q S dT - f Q S dT + [ Q S dT 

Jfc J ad Jad 

+ f Q dT + f QdT 
Jco Joa 


- [ Q S dT + [ Q S dT + [ Q S dT 

Jdefcbad Jfc Jad 

+ f QdT + [ Q dT 
Jco Joa 


J Xk — 


I DEFCBAD 


Q S dT + ( Jx k )n 


In Eq. (9), the first term on the right hand side is an integral on the closed 
contour DEFCBAD that excludes the crack tip. The second term on the right 
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hand side of the Eq. (9) corresponds to the line integrals on the crack faces. As 
will be shown below, the closed contour integral can be converted to an area or a 
domain integral. Hence, J Xk can be written as 

*^ x * — domain (*^ z *)/»ne (10) 

Closed Contour Integral — Invoking the divergence theorem, the closed con- 
tour integral of Eq. (9) can be converted to a domain integral as 


-L 


DEFCBAD 


q s d r = 


-L 


a 

[Wn k - Oij — ^ rij } SdT 
DEFCBA ox J 


= -/ 


f d(W5) d{oi JA 

[— o 5 dA 

A dx k dxi 


( 11 ) 


Hence, the domain integral is 


(7 \ - _ ( \ u/ dS du i dS ^ j a 

( **) domain J\ W ^ ^7 1 dA 


f . dW dta , 

A 1 a** " °'>dT k 1 s iA 


( 12 ) 


In deriving Eq. (12) the equations of equilibrium 



and the strain-displacement relationships 


e ij 


1 dui du J_ 

2 ' dxj dx{ 


were used. 
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In conventional finite element analysis, the equations of equilibrium are not sat- 
isfied pointwise in the domain that is modeled. Numerical experimentation showed 
that the differences between including and not including the terms involving the 
equations of equilibrium are of the order of 10 -3 to 10 -4 of the integral values for 
several problems. Therefore, in writing Eq. (12) the equations of equilibrium are 
assumed to be satisfied exactly. 

The term in brackets in the second integral in Eq. (12) is identically equal to 
zero pointwise for a linear elastic material. This term, however, is non zero in 
elastic-plastic problems. Although only linear elastic materials are considered in 
this report, this term is included because the algorithm is intended for use with 
nonlinear material problems. 

The domain integral in Eq. (12) can be rewritten in a form convenient to 
numerical computation as 


{J% k ) 


domain 


~L l 


w jr k - K*> r kK s '}] dA 


- J/£- {a}T{e '’> }]s iA 


(13) 


where 


w 


T 


{ <*11 <712 «22 } 



/ 

1 dx k 


2^£l2 d<22 •. 

dzfc dxk 



i 

1 dx k 


du2 

dx k 


} 


a 


<7n <712 

<712 <722 


{S'} T = { 


dS 

dx\ 


dS 
dx 2 ' 


(14) 
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The strain-energy density W is 


w — g ( CTllCl1 + ^ 22^22 + 2012*12] ( 15 ) 

The numerical implementation of Eq. (13) in a finite-element anlysis with 
isoparametric elements is presented in Appendix A. 

The integral J Xl in Eq. (13) for the linear elastic case is equivalent to the strain- 
energy release rate calculated by the virtual crack extension method since the right 
hand side of this equation represents the energy change per unit crack extension 
[8]. Equation (13) is also equivalent to an equation obtained by de Lorenzi [6,7], 
by calculating the change in energy by mapping a crack configuration with a crack 
length of a to another crack configuration of length a + da. 


Line Integrals — The four line integrals of Eq. (8) are 


(. Jx k )iine = [ QSdT + f QSdT + f QdT + [ Q dT (16) 

Jfc Jad Jco Joa 

When the term Q (defined by Eq. (2)) is zero on the crack faces, obviously, 
the line integrals are zero. On the crack faces, »i is always zero; when the crack 
face is traction free, a 12 and 022 are zero. Thus, for traction-free crack faces, the 
line integrals in Jxi always vanish. 

This is not the case with J X2 > the product integral. The term Q contains the 
term Wn 2 ; «2 = —1 on the line FO and ri 2 = 1 on the line OD. Thus, even 
for traction-free crack faces, Q (and line integral in J X2 ) can be non-zero. If only 
mode I deformations exist, then the strain-energy density W is zero on lines FO 
and OD, and, thus, Q and the J X2 line integrals are zero for traction-free crack 
faces. When mode II deformations occur and only the square-root singular term 
exists in the stress field, the sum of the integrals on FO and OD is zero due to 
the antisymmetric nature of the deformations, and again the J X2 line integrals are 
zero. But in a general crack problem subjected to external loading, both singular 
and non-singular stress fields exist around the crack tip; thus, the J %2 line integrals 
will be non-zero (see Appendix B for details) . Evaluation of the line J X2 integrals is 
complicated by the singular stress field and the strain-energy density at the crack 


9 


tip because of the line integrals on CO and OA. Thus, a method of calculating 
the mode I and mode II components in mixed mode problems, the decomposition 
method, in which the J X2 line integrals vanish is presented in the next section. 


Separation of Modes in Mixed Mode Problems 


Direct method — The two integrals J Xl and J X2 of Eq. (10) can be used to 
establish the individual modes. If Jj and Jjj correspond to the pure mode I and 
pure mode II deformations, then [3,10] 

J. ti = Ji + Jll 

(17) 

Jx 2 = ~ 2 V Jl J II 

From Eq. (17), J/ and Jjj can be computed in terms of J Xl and J X2 as 


J/ = - [ %/ Jxi — J%2 + 

4 


J II — ^ [ \J *^*1 ~ ^*2 “ 

J total = Ji + Jll 

For general mixed mode deformation the computation of J Xl and J X2 from 
Eq. (13) and the use of Eqs. (17) and (18) completely define the individual modes. 
However, for either pure mode I or pure mode II deformations the values of J Xl 
and J x 2 alone are insufficient to determine the individual modes. (Additional 
information such as the local crack tip deformation is needed to define the mode.) 
This is because in either case J X2 = 0 and J Xl gives the total integral. Equation 
(15) then suggests that only Jj exists. This difficulty is due to the choice of the 
positive sign for the square root terms in Eq. (18). 

Decomposition method — As previously mentioned, the advantage of trans- 
forming the contour integral into a domain integral is lost because of the non-zero 
line integrals in J X2 . These line integrals are necessary to account for the terms 
containing the product of the singular and non-singular stress (strain) fields in 


%/ Jx i + *f*2 
v^xTTJxT ] 2 
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the strain-energy density expression. It is shown in appendix B that the product 
terms can be eliminated by decomposing the stress and displacement fields into 
symmetric and antisymmetric parts. The resulting equation contains only the 
domain integral. This method termed as the decomposition method is described 
below. 

Consider two points P{x\,X 2 ) and P'(x 1 ,— £ 2 ) that tire in the immediate 
neighborhood of the crack tip and are symmetric about the crack line as shown 
in Figure 2. For general mixed mode deformations the displacements at P and P' 
can be expressed as a combination of symmetric and antisymmetric components 
as shown in Figure 2. Then 


and 


(“>'’} = (“ ls l + (19) 

l <*2P ) \ “2 S J l U 2AS ) 

1 u 1P , \ r «i S 1 r -« MS 1 (20) 

l “2 P' I ( -“2 S ) { "2AS ) 


where subscripts S and AS denote the symmetric and antisymmetric components, 
respectively. 

Equations (19) and (20) cam be used to determine the symmetric and antisym- 
metric displacements in terms of the displacements at points P and P' (see Figure 
2) as 



1 

2 

1 

2 


r u ip + u ip> 

\ U 2P ~~ u 2 P' 

f «ip - «ip' 
\ « 2 P + U 2P' 


} 

} 


( 21 ) 


Similarly, the symmetric and antisymmetric components of the stresses can be 
expressed in terms of the stresses at points P and P' (see Figure 3) as 
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> 



<711 1 

1 

<7llp - <7ll p , 1 

<722 | 

= - 
2 

| <722p - <722 p / \ 

<712 J 

AS 

[<712 P + <7l2p/ J 


( 22 ) 


The symmetric and antisymmetric displacements (Eq. (21)) and stresses 
(Eq. (22)) can be used to evaluate the four integrals JsxiiJsx 2 JaSxi, an <l JaSx 2 
using Eq. (10). Note that the integrals Js X2 an( i JaSx 2 (domain and line compo- 
nents individually) will be identically zero because of the symmetric and antisym- 
metric nature of the stress and displacement fields. 

The individual modes, Jj and J/j, are now 


Jl 

= J Sxi 


J II 

= J ASxi 

(23) 

Jtotal 

= JSx i + JASxi 



Obviously, this procedure involves an additional step to evaluate the symmetric 
and antisymmetric components from Eqs. (21) and (22). However, the individual 
modes are directly available from the domain integrals Jsn and JaSx i- For pure 
mode I problems JaSxi = 0 an( l Jtotal — Jsx i- For pure mode II problems 
Jsx i = 0 an< l Jtotal ~ JaSx i* 


Crack-Face Pressure Loading 


In the above discussion, the crack faces were assumed to be stress free. When 
the crack faces are subjected to applied loading, additional terms need to be 
included in the domain integral formulation. 

Again consider equation (10) 
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{Jx k ) 

domain + (^k)line 


( 10 ) 


Jx k ~ 


When the crack faces axe loaded, the line integrals are (see Fig. 1(b)) 


M = f [Wn k + c n ^- + <7 22 §^| S ixt 

J F dx k dx k 

+ j D [l Vn t - <^|^] S dx 

+ / c > Wn * + + 

/O r„, ^Ul 5U2, , 

+ J a [ Wn k ~<712 q^ ~ ^22^^] 


(24) 


In deriving Eq. (24), ni = 0 and «2 = 1 were used on line AD and n\ = 0 and 
i %2 = — 1 were used on line FC. Note that when the crack faces are stress free 
<7 12 = <722 = 0. Therefore, all the line integrals vanish yielding only the domain 
integral of Eq. (10) as J Xk . 

For the decomposition method, the symmetric and antisymmetric line integrals 
Jsx k and J ASx k can be obtained by substituting the symmetric and antisymmetric 
components of stresses and displacements into Eq. (24). This yields 


(Jsxi)line 




dx i 


(JSxi)line — 0 


{JaSxi )line 




(25) 


{.J AS inline — 0 


Case 1 — When the only loading on the crack faces is a uniform pressure of 
magnitude p, 
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<712 = 0 ; °22 - -p 

The symmetry of the deformation requires that 


and 


du 2 1 _ j 

dx\ DA dx\ FC 


[£i 

dx 2 DA 


= ifr> 

OX 2 JTC' 


Therefore, the line integrals reduce to 


(26) 


(27) 


( Jxi)line — {Jsxi)line — 2 J ^ ^ Q X \ ^ ^ i l ^ dx\ ^ E 1 


(*^*2 )/me — (*^Sij)/me — {^ASxi)line ~ ^ASij)/*ne 0 


(28) 


where (~h,0) and (— /i,0) are the coordinates of points F and C in Fig. 1(b). 
Case 2 — When the only loading on crack faces is a uniform shear of magnitude 




<712 = T 5 <722 = 0 


(29) 


The antisymmetry of deformation requires that 


[£i = -i£i 

<7*1 ZM <7^1 FC 


and 
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Therefore, the line integrals reduce to 


(*^*i)fme — {^ASxi)line ~ 2 J 

(*^*3 )ltne = ( )itne = (^Sxi 


RESULTS AND DISCUSSION 

In this section, the EDI method is applied to two mode I problems and two 
mixed mode problems. Because only linear elastic materials are considered here, 
to validate the EDI method the J— domain integrals are compared to the strain- 
energy release rate values from the literature [14-18]. Unless otherwise specified, 
a Poisson’s ratio of 0.3 was used for all configurations analyzed. 

Mode I problems 

Remote loading — The first two examples considered are the center cracked 
tension (CCT) and single-edge cracked tension (SECT) specimens as shown in 
Figure 4(a) and (b), respectively. A crack-length-to-width ratio of 0.8 was 
considered. Due to the symmetries in the problems, either one-quarter or one- 
half of the specimen was modeled with 34 8-noded parabolic elements as shown 
in Figures 4(c) and (d). At the crack tip, collapsed quarter-point elements were 
used. Figure 5 shows the detail near the crack tip and also the domains D\,D 2 , 
and Z >3 used to evaluate the domain integrals. 

Figure 5 presents a variety of S-functions that were used in computing the 
domain integrals. Type I is a simple linear function. Types II and III are quadratic 
functions. These S-functions were created by defining the values of S at nodes on 
the elements, as discussed in appendix A. The functional forms corresponding 
to the S-functions are given in Figure A2. These three functions were used in 
conjunction with each of the domains £>i, £> 2 , and £>3. The other two S-functions 
in Figure 5 involved two layers of elements, i.e., domains D\ and £> 2 , or £>2 and 


du\ f° du\ 

r_i s dxi + 2 r— — dx 1 

dxi J_i dxi 


(31) 


£>3* These two functions were used to study the influence of using more than one 
layer of elements in the computation of the domain integrals. 

Table 1 presents the normalized values of J Xl (Jsxi an< l J i) obtained by 
considering various domains and S-functions for both CCT and SECT specimens. 
(Due to symmetry, the J Xl values in this table correspond to twice the value 
computed with the elements in the various domains shown in Fig. 5. Also note 
that the J Xi values computed from the elements in the domains D{ shown in Figure 
5 were non-zero. However, the symmetric group of elements below the crack line 
contribute a value with the same magnitude but with opposite sign. Hence, the 
total value of J X2 is equal to zero, as required by these mode I problems.) For 
comparison the normalized values of the strain-energy release rate G (and hence 
J Xl or Jsxi) from the literature [15-18] and obtained by the virtual crack closure 
technique (VCCT) [14] using the same finite element model are included in this 
table. For the CCT specimen, the various S-functions and domains give nearly 
identical results; the variation of the computed values was about 1.2 percent. For 
the SECT specimen, the domain integrals were nearly same with a maximum 
variation of 2 percent. The values computed using the EDI are, as expected, in 
excellent agreement with the VCCT calculations. The results in this table suggest 
that domains consisting of a single layer of elements and linear S-functions are 
sufficient to obtain accurate results. 

The model in Figure 4(d) is locally symmetric about the crack tip. As 
mentioned previously, the finite element model need not be symmetric about a 
line normal to the crack line and parallel to the a^-axis. To verify this, a new 
model without symmetry about the crack tip was developed. A simple rectangular 
idealization was used throughout the model as shown in Figure 6. Note that in this 
model singularity elements were not used at the crack tip. Three domains were 
used in the EDI calculations as shown in Figure 6. Linear S-functions were used 
in the EDI calculations. Table 2 presents the normalized values of J Xl obtained 
with the three domains. All domains gave nearly identical results; the variation 
of the computed values was about 0.6 percent. Furthermore, these results agree 
extremely well with those calculated with quarter-point singularity elements at 
the crack tip and with the VCCT method. The normalized value of the integral 
with this model without singularity elements is about 3 percent lower than the 
value with singularity elements. Note that the application of the VCCT method 
is difficult when the finite-element model is not locally symmetric about the crack 
tip. The EDI method does not have such a restriction. 

Crack-face loading — As pointed out earlier, the equivalent domain integrals 
need to be modified when the crack faces are subjected to external loading. 
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Consider the CCT specimen in Figure 4(a). Instead of the remote tensile loading 
shown, let the specimen be subjected a uniform crack-face pressure of magnitude 
p. Figure 7 presents the domains near the crack tip for one or two layers of 
elements. For each of the S-functions the line integrals in Eq. (28) were computed 
and are shown in Table 3. 

Table 4 presents the domain and line components of the integrals J Xl ( J Xl — 
Jsxi = Jl) f° r the crack-face pressure loading. The sum of the domain and line 
components computed for the crack-face pressure loading is exactly equal to the 
domain integrals computed for the remote loading situation. This is true for all 
domains. 


Mixed Mode Problems 

The second two examples considered are mixed mode problems, an angle crack 
problem (Fig. 8(a)) and a cracked bimaterial problem (Fig. 8(b)). These problems 
are used to discuss the separation of individual modes using the domain integrals. 

Angle crack — Figure 8(a) shows the configuration used for an angle crack in 
a finite plate. Figure 9(a) and (b) show the finite element idealization with 77 
8-noded parabolic elements. At the crack tip, collapsed quarter-point elements 
were used. 

Four domains D\ and Z?2 (Figure 9(b)) were used to evaluate the domain 
integrals. Three S-functions, Types I, II, and III, were used with the two domains. 
As mentioned previously, Jj and J/j can be evaluated by (a) computing J Xl and 
J X2 and using Eq. (18) or by (b) computing Jsx k an d / ASx k and using Eq. (23). 
Table 5 presents the results obtained using both of these methods. 

First, note that the total integral J Xi is identically equal to the sum of J Sx i 
and J ASx\' Next, the integrals Js Xl («//) and JaSxi i^Il) agree very well with 
the Jj and Jjj from the literature [14-18], Third, the two domains and three 
S-functions showed a maximum variation of 0.3 percent for Js Xl and 1.5 percent 
for JaSxi • Thus J Xl from the undecomposed displacement field and Js Xi aud 
JaSxi from the symmetric and antisymmetric displacement fields yield accurate 
results. 

The J// values calculated using J Xl and J X2 are not accurate. The errors are in 
the calculation of the J X2 integral. As previously mentioned, J X2 includes the line 
integral in the singularity region and hence requires integration of singular terms 
near the crack tip. The line integrals in the singularity element region (i.e. in 
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quarter-point element region) were computed in the same manner as for non- 
singular elements, using a 2-point Gaussian quadrature. This may be the reason 
for about 3 percent error in the J — values in Table 5. When J X2 was evaluated 
without including the line integrals very large errors were found and the results 
were domain dependent. Therefore, the line integrals are necessary for accurate 
evaluation of J X2 . 

The integrals Js Xj and JaSx 2 are also product integrals like J X2 . These 
integrals were very nearly zero (within machine accuracy) for all domains and 
for all S-functions. Note that Jsx 2 [JaSx 2 ) was evaluated with the symmetric 
(antisymmetric) components of the displacements and stresses. Hence, the zero 
values computed for Jsx 2 and JaSx 2 were expected because of the symmetric and 
antisymmetric nature of the displacement field used in the computations. Also 
note that the line integral components are identically zero. From these results 
the decomposition method appears to be the preferable over the direct method 
(Eq. 18) to separate the modes. 

Cracked bimaterial plate — The second mixed mode problem analyzed was 
a bimaterial plate with central crack along the interface (Fig. 8(b)). The plate 
was subjected to remote tensile loading. Different combinations of materials were 
used in the plate. The materials and their properties are summarized in Table 6. 
Because of symmetry, only one-half of the plate was idealized as shown in Fig. 10(a). 
Figure 10(b) shows the detail near the crack tip for a model with ^=0.05. This 
model had 517 nodes and 156 elements. At the crack tip, quarter-point singularity 
elements were used. Four domains and linear S-functions were used in the domain 
integral calculations. 

In contrast to the angle crack problem where only the displacement field was 
decomposed into symmetric and antisymmetric components, in this problem both 
the displacements and the stress fields were decomposed. This was done because in 
the bimaterial problem, the stresses derived from a symmetric (or antisymmetric) 
displacement field yield both symmetric and antisymmetric (or antisymmetric and 
symmetric) stresses. The decomposed stresses and displacements were used to 
calculate the domain integrals. 

Table 7 presents the integrals Jg Xl and JaSxi f° r various combinations of 
materials and domains considered. These integrals are also compared to the strain- 
energy release rate calculated using VCCT formulae given in reference 13. First, 
note that the total integrals show less than one percent variation among the four 
domains for all material combinations. The total integrals also agree extremely 
well (within 1.6 percent) with the total strain-energy release rate calculated with 
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VCCT. Next, the mode I integrals show some variation among the four domains. 
The maximum variation between D\ and £>4 is 4 percent for the aluminum-epoxy 
case, and is 4.8 percent for the steel-epoxy case. The differences between the 
mode I integral values of domain D\ and the VCCT values are 3, 3.4, and 0.2 
percent for the aluminum-epoxy, steel-epoxy, and steel-aluminum combinations, 
respectively. The discrepancies are much smaller with other domains. In contrast 
to the total and the mode 7, the mode 77 integrals show large variations among all 
the domains for all material combinations. (The variation of the mode 77-to-total 
ratio is about 5 percent between the domains D\ and £>4 for the steel-epoxy case.) 
Furthermore, the mode II integrals do not agree with those calculated with the 
VCCT. Also note that the discrepancies are much larger for the steel-epoxy case 
than the other material combinations. 

The mode 7 and mode II domain dependence may be related to the oscillatory 
singularity that exists at the crack tip in a bimaterial plate [19]. The VCCT also 
yields inconsistent values of individual modes for different values of A, the crack 
closure length [19]. To determine the effect of A, the crack-tip singularity element 
size was reduced from ^ = 0.05 to 0.025. The new model with ^ =0.025 is shown 
in Fig. 10(c). This model was obtained by subdividing each of the 8 singularity 
elements into two elements; a regular 8-noded element and a singularity element. 
This new model had 543 nodes and 164 elements. Five domains were used with 
this new model but the domains D\ to D 4 were the same as the model of Fig 10(b) . 
Because the steel-epoxy has the largest differences in the material properties (and 
hence the largest oscillatory power), the steel-epoxy case was analyzed with this 
new model. 

Table 8 presents the mode I, mode II, and total integrals obtained with various 
domains and with the two models of Figs. 10(b) and 10(c). The results from the 
two models are almost identical for domains D\ to £>4. The total integral for 
domain D 5 is only different by 0.3 percent when compared to that of domain £>4. 
Thus, the model with ^ =0.025 shows the same domain dependence of the mode 
I and mode II values. The mode I and mode II VCCT results, as expected, also 
show a dependence on A, the crack closure length. In contrast, the total strain- 
energy release rate is virtually unchanged as £ is reduced from 0.05 to 0.025. As 
shown in reference 19, the mode I and mode II strain-energy release rates do not 
have well defined limits as A approaches zero. However, the total strain-energy 
release rate is independent of A. This is exactly the behavior displayed here by 
the J— integrals. Therefore, the inconsistent behavior of mode I and mode II 
integrals is due to the oscillatory singularity. The different domains correspond 
to different values of A in the VCCT. Hence, the domains D 1, £>2, etc. show a 
continuously changing value of mode 7 and mode II while the total J remains 
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unchanged. Thus, the EDI method and VCCT are consistent and are simply 
reflecting the continuum behavior. 

CONCLUDING REMARKS 

An equivalent domain integral (EDI) method for the calculation of J — integral 
fracture parameters (total and product integrals) is presented for two-dimensional 
elastic cracked bodies. The integrals consist of an area or domain and line integrals 
on the crack faces. The line integrals vanish only when the crack faces are traction 
free and the loading is either pure mode I or pure mode II or a combination of 
both with only the square-root singular term in the stress field. However, for finite 
size crack bodies subjected to mixed mode loading, even when the crack faces are 
traction free, the line integrals on the crack faces are non zero and are necessary 
to account for the non-singular field. Two methods for calculating the J— domain 
integrals for individual modes are discussed. 

The EDI method gave accurate values of the J— integrals for the various mode 
I and mixed mode problems studied. The calculated integrals agreed very well 
with the strain-energy release rates from the literature. The numerical studies 
showed that domains consisting of one layer of finite elements were sufficient to 
give accurate results. Several types of S-functions, the functions that facilitate 
the conversion of the line J— integrals to the domain integrals, were used in the 
numerical studies. These studies showed that the results are independent of the 
type of S-functions studied and even the use of simple linear S-functions yields 
accurate results. 

Two methods, direct and decomposition, for evaluating the individual mode 
components were studied. In the direct method the total and product J— integrals 
are computed and the individual modes are calculated from these two integral 
values. In the decomposition method, the displacements and stresses are first 
decomposed into the symmetric and antisymmetric components and then the 
J— integrals are evaluated for each of these components. The decomposition 
method was found to be more accurate than the decomposition method. The 
direct method requires evaluation of line integrals on the crack faces and hence 
requires integration of singular strain energy density near the crack tip. Numerical 
integration of these terms contributed to the loss of accuracy in this method. In 
contrast the line components are identically zero in the decomposition method. 

For cracks in a homogeneous material only the displacement field needs to be 
decomposed into the symmetric and antisymmetric components. On the other 
hand, for bimaterial plates with cracks along the interface both the displacement 
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and stress fields need to be decomposed. For a bimaterial specimen the domain 
integrals calculated with the decomposed displacement and stress fields produce 
mode I and mode II J — integral values that show domain dependence in a 
manner similar to the dependence on the crack closure length in the virtual crack 
closure technique (VCCT). This behavior is caused by the oscillatory part of the 
singularity in the bimaterial crack problem. However, the total J— integral is 
domain independent. 
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APPENDIX A 

NUMERICAL IMPLEMENTATION OF THE EDI METHOD 


This appendix presents the numerical implementation of the EDI method in 
the context of a finite element analysis with isoparametric elements. Some of the 
details presented in this appendix can also be found in references 8,10, and 11. 

A typical finite element model with 8-noded isoparametric elements that is 
locally symmetric about the crack line (x 2 = 0 ) is shown in Figure Al. This 
type of locally symmetric model is necessary for the computation of Jsx% > ^ASxi i 
etc. (Eq. (19-22)). In contrast, either symmetric or non-symmetric models can be 
used if J Xl and J X2 are to be computed. Note that symmetry of the model with 
respect to the x\ = 0 line is not a requirement. 

Consider the two contours To and Ti that constitute one layer of elements as 
shown in Figure Al. More than one layer of elements may be considered. However, 
one layer is considered here for convenience in presentation. As is shown in the 
results and discussion section, the results do not depend on the number of layers 
of elements between To and Ti and one layer of elements is computationally very 
convenient. 

The elements I, J,K, L,L / ,K < , J*, and I* are bound by the contours To and Tj. 
Thus the domain integral is 


I,J, • i' 

(*^ x k) domain ^ v (-^1) 

t 

where J Xk , is the area integral over the i th element. 

In the isoparametric representation, the displacements within the element are 
defined by the shape functions Nj and the nodal displacements (u a )j 


u a = Nj{u a )j 


(A2) 


where a =1,2, Nj = Nj(^,rj), £,r} are the parent coordinates and j = 1 to the 
number of nodes on the element. 
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The area integral J Xk . of the i th element, Eq. (13), can be computed using 
Gaussian quadrature as 


- - [ LL [w £i - { “*‘ }T |2|<S ' } 1 {det |j|) « M 


+ J J i - {*} T (4* > 1 s ( det [ J ]) d z dT > ]< 


Mg 

= - E 


m=l 


M g 

E 

n=l 




+ [^— — {o } T {ei fc }! s } u) m UJrvfdc([J]) | ^ 


(A3) 


where Mq is the number of Gaussian quadrature points used in each direction £ 
and T), w m and w n are the Gaussian weights and det [J] is the determinant of the 
Jacobian matrix [J] defined by 


[J] 




'(ll)l (X2)l‘ 

• ... dNs - 


(ll)2 (X2)2 

dNj dNi dNs 

. dtf drf drf . 


-(®l)8 (^ 2 ) 8 - 


(A4) 


Most of the quantities necessary for Eqs. (13-15), W , {a}, {ei fc }> {u* fc }> an ^ [fi]* 
are readily known in terms of the nodal displacements of the element. However, 
computation of the terms S, {S'}, and needs special attention and is discussed 
below. 

S-Functions — As mentioned previously, the S-function needs to be continuous 
and have values equal to zero on the contour T i and unity on the contour Tq. The 
S-function can be conveniently defined using the element shape functions as 
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S(C,r,) = Nj Sj 


(/ 15 ) 

where j = 1 to the number of nodes on the element and Sj is the nodal value 
of the S-function at node j on the element. The requirement that the function 
be zero on Ti and unity on To can be satisfied by specifying zero values for S at 
nodes 1, 8, and 7 on the rj — 1 line and unit values at nodes 3 , 4 , and 5 on the 
r] = — 1 line of the parent element (see Figure Al). All values of Sj except those 
at the midside nodes 2 and 6 are now defined. Various choices of Sj for nodes 2 
and 6 lead to different types of S-functions and these are illustrated in Figure A2. 
Thus specifying the nodal values Sj completely defines the S-functions. In fact a 
single subprogram that defines the nodal values Sj and, hence, the S-function can 
be used for all elements between Tq and Ti. The only requirement for all elements 
I, J, • • • l' is that the lines with 77 = —1 should correspond to the contour To 
and the lines with r? = 1 should correspond to the contour T \ ( see elements K 

f ' 

and K in Fig. A2). This can be achieved by defining the element connectivity 
appropriately. 

q C 

Partial Derivatives of S — Once S is defined the partial derivatives of S, 
and can be easily computed using the isoparametric formulation as 




(A6) 


where [J] is the Jacobian matrix defined in Eq. (A4). 

Partial Derivatives of W — The terms are computed using a bilinear 
extrapolation (in terms of the parent coordinates f and rj) for W and using the 
values of W at the 2x2 integration points. In reference 10, the integral f fj-dA 

nri7 

was approximated with a one-point integration by evaluating at the center of 
the element. A slightly different approach is taken here. Because all the quantities 
are known at the integration points, the integration is carried out in the same 
manner, without further approximations, as for the other terms in Eq. ( 13 ). The 
values of the stresses are known to be more reliable at the 2x2 Gaussian points 
within the element (in comparison to the nodal values). The strain- energy density 
W is approximated in the bilinear form as 

W[£,ri) =01+ (A 7 ) 
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Using the 2x2 Gaussian values of the strain-energy density W, Eq. ( A7 ) can be 
rewritten as 


where 


and 


Wte>v) =[!*»» Cti][T}{W g } 


1 

1 

-Vs 

1 

-vs 

1 

V3 

1 

vs 

4 

-Vs 

_ 3 

Vs 

-3 

-vs 

-3 

Vs 

3 


{W g } t ={W> W n w m Wrv } T 


(AS) 


(A 9) 


(-A10) 


Wl,Wu,Wui t and Wjy are the values of W at the 2x2 Gaussian points shown 
in Figure A3. The partial derivatives and are 


( aw >, 
I -5T 1 


'o l o V 

1 aw 1 
l -&o ) 


° 0 1 


[T] [W a \ 


(^ 11 ) 


Equation (All) simplifies considerably and can be rewritten as 


dW 

dZ 



(W T + W n - W UI -Wjy ) 


+ -r (Wi - w n - w HI +wjv ) 

4 


dW 
d V 


~ (Wi - Wn + Wm -Wrv ) + ^ ^ ~ Wn ~ Wfn +W/V ) 

(A12) 


The derivatives can now be obtained as 
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/ dW \ 

1 

aw | 

V 3X2 / 




(X13) 


where [J] is defined in Eq. (A4). 

All the necessary quantities in Eq. (13) are now known and, hence, the domain 
integrals can be calculated. 
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APPENDIX B 

NONSINGULAR STRESS FIELD IN THE COMPUTATION 

OF (Jx 3 )line INTEGRALS 

This appendix examins the ( J I2 )/, ne integrals when both the singular and 
nonsingular terms exist in the stress field. This appendix shows that line integrals 
vanish when the decomposition method is used. 

For the purpose of this discussion consider a 2-D cracked body subjected to a 
mixed-mode stress field ( Kj and Kjj) along with a uniform stress, o\\ = oq. The 
stresses along the two crack faces FO and DO (see Figure 1(b)) can be written in 
terms of the polar coordinates (r, 0) with the origin at the crack tip as [15] 


—2 Kjj 

*ll = + oo (Bl) 

\/2 ir r 

on FO and 


*n = + <r 0 (B 2) 

y/2nr 

on DO. The stresses 022 and <712 are identically zero on both FO and DO since the 
crack faces are stress free. Also note that the mode I stress field does not contribute 
to the 0 n stresses on either FO or DO. Assuming plane strain conditions, strains 
eil are 


on FO and 


en 


E L y/2 IT r 



€ 11 


(l-i/ 2 ) r 2 
E \-y/2irr 



(BZ) 


(B 4) 


on DO. 

For traction-free crack face problems Q, see equation (2), reduces to 
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(B 5 ) 


Q = W n k = - an e n n k 

where W is the strain energy density and n k is the component of the vector normal 
to the crack face in the k tfl direction. 

The line integrals corresponding to J Xl vanish because n\ = 0 on FO and OD 
and hence Q is identically zero. Therefore, only the (*/x 2 )/, nc integral needs to be 
examined. Using Eqs. (Bl) - (B5) and assuming that CO = OA, and FC = DA, 
the ( J X2 )j t - ne integral can be written as 


(1 - „*) 



2 Ku )2 
s/YtTt 

\2 

s/2 ~1 Tt 


S dr -(- 


S dr + 


L 

Fo 


r D 

ta 

r A 


{°o + 
(^0 + 


2 Kn 
V2 tTt 


) 2 S dr 


2 Kji 
y/2 7 r r 


) 2 S dr 
(B 6) 


where rp = rjj = /2 and rc = r a = h- Note that (— /2>0) an< i ( — /i, 0) are the 
coordinates of the points F and C in Figure 1(b). The line integral of Eq. (B 6) 
simplifies to 


(Jx 2 h 


4(1 -u 2 ) 


me 


E 


<TQ Kjj 


r f rF s 

f rc 1 1 

/ dr + 

/ dr 

\-Jrc v2 7r r 

lo \J2-Kr J 


( 57 ) 


The contributions in Eq. (57) are due to the product of singular {Kjj) and 
nonsingular ((To) terms. If only a pure singular stress field exists in the cracked 
body (i.e., when oo = 0), then the (•fx 2 )/ tne is identically equal to zero. Therefore, 
the presence of a nonsingular stress field in a mixed-mode problem gives a non 
zero (*/x 2 )/ tne integral. Numerical evaluation of this integral is difficult near the 
crack tip due to 1/y/r singularity (see Eq. 57). 

Using the decomposition method, the symmetric and antisymmetric stress 
fields along the crack faces can be written as 

Symmetric : 
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(^llJsinyu/ar = 0 

(*^11 ) Nonsingular ^0 


(Bll) 


Antisymmetric : 


(^11 ) Singular 


— 2 Kn 

y/2lrr 


(^ll ) N onsingtdar ® 


(^12) 


Note that the product of singular and nonsingular terms for both symmetric 
and antisymmetric stress fields is identically zero. From Eq. (B 7), it is obvious 
that the line integrals for conditions ( B8 ) and (59) vanish separately. Because 
the decomposition method considers symmetric and antisymmetric components 
separately, their cross terms are eliminated. Thus, ( Jsx 2 an< l JaSx 2 ) line integrals 
are identically equal to zero. 
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Table 1 Comparison of non-dimensional J-integral for CCT and SECT speci- 
mens subjected to remote tensile loading. (Plane strain, j = 0.8) 


Domain 

S-Function 

EJ XI 

p 2 7ra(l — 

CCT 

SECT 

D i 


3.285 

137.1 

D 2 

Type I 

3.284 

136.7 

D 3 


3.283 

137.1 

Di 


3.274 

135.8 

D 2 

Type II 

3.273 

136.4 

D 3 


3.270 

137.0 

£>i 


3.297 

138.4 

^2 

Type III 

3.295 

137.0 

D 3 


3.296 

137.2 

DxkDi 

Type IV 

3.302 

137.0 



3.309 

136.6 


Type V 

3.290 

135.7 

D2&D3 


3.298 

136.3 


VCCT [14] 

3.396 

137.4 


Reference Value [15-17] 

3.298 

143.8 
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Table 2 Comparison of non-dimensional J-integral for CCT specimen sub- 
jected to remote tensile loading computed with a model with non- 
singular elements. (Figure 6, Plane strain, y = 0.8) 


Domain 

EJ XI 

p 3 7ra(l— v 2 ) 

D i 

3.201 

D 2 

3.221 

D 3 

3.219 
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Table 3 Line integrals for various S-functions (see Figure 7) 


One layer of elements in the domain 

S- Function 

Line-integrals (J x \)u nt and (Jszi)iine 

Type I 

= l,Sj =0.5,5* = 0) 

To [ u 2,- + 4 u 2> + u 2fc ] 

Type II 

(Si = 1 ,Sj = S k = 0) 

^[3u 2< + 4u 2 . - u 2k ] 

Type III 

(St = Sj = lySk — 0) 

i + 4u 2j . + 3u 2fc ] 

Two layers of elements in the domain 

Type IV 

(Si = Sj = S k = 1, 5, = 0.5, S m = 0) 

^[u 2fc + 4 u 2 , + u 2m ] 

Type V 

(S i = S j = S k = l,S l = S m =0) 

^[3u 2fc + 4 u 2 , - u 2m j 
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Table 4 Comparison of non-dimensional J-integral for crack face pressure and 
remote loadings for the CCT specimen. (Plane strain, £ = 0.8) 






EJ Xi 

p a 7T a(l — 


Domain 

S- Function 

Crack face pressure loading 

i.^xi)domatn i^xi)line xi)total 

Remote loading 

(Jx, ) domain 

D i 


1.872 

1.413 

3.285 

3.285 

D 2 

Type I 

2.064 

1.220 

3.284 

3.284 

D 3 


2.322 

0.961 

3.283 

3.283 

D i 


1.916 

1.357 

3.273 

3.273 

D 2 

Type II 

2.127 

1.146 

3.273 

3.273 

D 3 


2.412 

0.858 

3.270 

3.270 



1.828 

1.469 

3.297 

3.297 

£>2 

Type III 

2.001 

1.294 

3.295 

3.295 

Dz 


2.233 

1.063 

3.296 

3.296 

D\iiDz 

Type IV 

1.889 

1.413 

3.302 

3.302 

Dz&zDz 


2.089 

1.220 

3.309 

3.309 

D 1 &D 2 

Type V 

1.933 

1.357 

3.290 

3.290 

Dz&tDz 


2.152 

1.146 

3.298 

3.298 
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Table 5 Comparison of frg, ra {f_ t ,gyj computed from direct and decomposition 
methods for an angle crack problem. (Figure 7, Plane strain,^ = 
0.5, <j> = 45°) 






Direct method 

Decomposition Meth 

Domain 

S-Function 

Jxi 

Jx 2 

h 

J II 

•W Ji) 

JaSxi ( J II ) 

D i 


1.781 

-1.317 

1.489 

0.291 

1.439 

0.342 

D 2 

Type I 

1.774 

-1.304 

1.481 

0.293 

1.437 

0.338 

D 3 


1.775 

-1.320 

1.481 

0.294 

1.437 

0.338 

Da 


1.775 

-1.324 

1.478 

0.296 

1.437 

0.338 

D i 


1.781 

-1.320 

1.488 

0.293 

1.440 

0.341 

D 2 

Type II 

1.774 

-1.318 

1.480 

0.294 

1.436 

0.338 

Dz 


1.774 

-1.320 

1.479 

0.294 

1.436 

0.338 

Da 


1.770 

-1.321 

1.474 

0.296 

1.432 

0.338 

D i 


1.780 

-1.314 

1.491 

0.290 

1.438 

0.343 

£>2 

Type III 

1.775 

-1.318 

1.482 

0.293 

1.437 

0.338 

Dz 


1.776 

-1.320 

1.482 

0.294 

1.438 

0.338 

Da 


1.779 

-1.326 

1.482 

0.297 

1.441 

0.338 





VCCT 

[14] 

1.447 

0.340 





Collocation 

[15-17] 


0.325 





BFM 

[18] 

1.416 

0.330 
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Table 6 Material properties used in the bimaterial plate problems. 


Material 

E, psi 

V 

Epoxy 

0.45 x 10 6 

0.35 

Aluminum 

10.0 x 10 6 

0.30 

Steel 

30.0 x 10 6 

0.30 
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Table 7 Comparison of the non-dimensional J-integral for a center cracked 
bi-material plate subjected to remote tensile loading. (Plane strain, 
1 = 0 . 1 ) 



Material 



EtJ 

P 3 ’ ro(l — j/j*) 



Domain 



1 

2 


Mode I 

Mode II Total 


Aluminum 

Aluminum 

D i 

1.012 

0.0 

1.012 



D 2 

1.012 

0.0 

1.012 



D 3 

1.012 

0.0 

1.012 



D 4 

1.012 

0.0 

1.012 




(1.010)* 

(0.0) 

(1.010) 


Aluminum 

Epoxy 

D i 

8.771 

0.590 

9.361 



D 2 

8.659 

0.716 

9.375 



D 3 

8.571 

0.807 

9.378- 



Da 

8.426 

0.989 

9.415 




(8.671) 

(0.613) 

(9.284) 


Steel 

Epoxy 

D i 

24.44 

1.924 

26.37 



d 2 

24.09 

2.329 

26.41 



D 3 

23.81 

2.618 

26.42 



D 4 

23.34 

3.204 

26.54 




(24.13) 

(1.996) 

(26.13) 


Steel 

Aluminum 

D i 

1.893 

0.0333 

1.931 



D 2 

1.885 

0.0467 

1.931 



Dz 

1.875 

0.0567 

1.932 



Da 

1.859 

0.0750 

1.934 




(1.863) 

(0.0581) 

(1.921) 
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Table 8 Comparison of the non-dimensional J-integral for a center cracked 
steel-epoxy bi-material plate subjected to remote tensile loading for 
two values of (A/a) = 0.05 and 0.025. (Figure 8, Plane strain, 
f=0.1) 


A 

Domain 


EiJ 

p a ira(l— v'{) 


a 


Mode I 

Mode II 

Total 



D i 

24.44 

1.924 

26.37 


D 2 

24.09 

2.329 

26.41 

0.05 

D 3 

23.81 

2.618 

26.42 


D 4 

23.34 

3.204 

26.54 


VCCT 

(24.13)* 

(1.996) 

(26.13) 



D i 

24.44 

1.924 

26.37 


D 2 

24.09 

2.325 

26.42 

0.025 

D 3 

23.81 

2.614 

26.42 


D 4 

23.31 

3.228 

26.53 


D 5 

22.62 

3.983 

26.61 


VCCT 

(23.43) 

(2.765) 

(26.19) 


* Values in the parentheses were obtained with VCCT formulae given in reference 14. 
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Figure 1 - Contours around the crack tip. 
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Figure 2 - Symmetric and antisymmetric 
displacement components. 



I 
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Figure 3 - Symmetric and antisymmetric stress 
components. 
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Figure 4 - Mode I problems and the finite element 
idealization. 
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Figure 5 - Various domains and S-functions used in mode I 
problems. 





(a) Domain with one layer of 
elements. 



(b) Domain with two layers of 
elements. 


Figure 7 - Nodes used in calculating the line-integrals 
for crack face pressure loading. 
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Figure 8 - Mixed mode problems. 
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Figure 10 - Finite element idealizations and domains used 
for the bimaterial plate problem. 
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Figure A1 - Typical finite element model. 



Parent Element S-Functions 
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Figure A2 - Various S-functions. 
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Figure A3 - 2x2 Gaussian quadrature points in 
the parabolic element. 
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16. Abstract 

An equivalent domain integral (EDI) method for calculating J-integrals for two-dimensional cracked elastic 
bodies is presented. The details of the method and its implementation are presented for isoparametric 
elements. The total and product integrals consist of the sum of an area of domain integral and line integral;; 
on the crack faces. The line integrals vanish only when the crack faces are traction free and the loading is 
either pure mode I or pure mode II or a combination of both with only the square-root singular term in the 
stress field. The EDI method gave accurate values of the J-integrals for two mode I and two mixed mode 
problems. Numerical studies showed that domains consisting of one layer of elements are sufficient to 
obtain accurate J-integral values. Two procedures for separating the individual modes from the domain 
integrals are presented. The procedure that uses the symmetric and antisymmetric components of the stres ; 
and displacement fields to calculate the individual modes gave accurate values of the integrals for all 
problems analyzed. The(LDljnethod when applied to a problem of an interface crack in two different 
materials showed that the mode I and mode 13 components are domain dependent while the total integral is 
not. This behavior is caused by the presence of the oscillatory part of the singularity in bimaterial crack 
problems. The EDI method, thus, shows behavior similar to the virtual crack closure method for bimateri: 1 
problems. 
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